Figure 4

Load libraries

library(STutility)
library(ggplot2)
library(ggpubr)
library(dplyr)
library(magrittr)
library(UpSetR)

Assemble spaceranger output files and merge curated meta data

# Mouse brain
samples <- Sys.glob(paths = "../data/spaceranger_output/smallintestine/*/filtered_feature_bc_matrix.h5")
imgs <- Sys.glob(paths = "../data/spaceranger_output/smallintestine/*/spatial/tissue_hires_image.png")
spotfiles <- Sys.glob(paths = "../data/spaceranger_output/smallintestine/*/spatial/tissue_positions_list.csv")
json <- Sys.glob(paths = "../data/spaceranger_output/smallintestine/*/spatial/scalefactors_json.json")

infoTable <- data.frame(samples, imgs, spotfiles, json)
infoTable <- cbind(infoTable, arrayID = do.call(rbind, strsplit(infoTable$samples, "/"))[, 5])

curated_metadata <- openxlsx::read.xlsx("../data/sheets/curated_RNA_rescue_sample_metadata.xlsx", sheet = 3)
curated_metadata <- setNames(curated_metadata, nm = c("storage_time", "seq_date", "include", "RNA_rescue",
                                                       "project", "experimenter", "processer", "comments",
                                                       "ID", "paper_id", "RIN", "DV200", "protocol", 
                                                      "source", "arrayID", "spots_under_tissue",
                                                       "genes_detected", "fraction_spots_under_tissue",
                                                       "median_genes_per_spot", "median_UMIs_per_spot", 
                                                       "saturation", "reads_mapped_to_probe_set",
                                                       "reads_mapped_confidently_to_probe_set",
                                                       "reads_mapped_confidently_to_filtered_probe_set",
                                                       "reads_mapped_to_genome",
                                                       "reads_mapped_confidently_to_genome",
                                                       "number_of_panel_genes"))

infoTable <- merge(infoTable, curated_metadata, by = "arrayID")

Load data into a Seurat object

SI <- InputFromTable(infoTable %>% arrange(seq_date))

## Using spotfiles to remove spots outside of tissue
## Loading ../data/spaceranger_output/smallintestine/V19T26-028_A1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V19T26-028_B1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V19T26-028_C1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V19T26-028_D1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-048_A1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-048_B1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-048_C1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-048_D1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-050_A1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-050_B1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-050_C1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10F24-050_D1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V10S29-108_A1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../data/spaceranger_output/smallintestine/V11B18-363_D1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## 
## ------------- Filtering (not including images based filtering) -------------- 
##   Spots removed:  21  
##   Genes removed:  6589  
## Saving capture area ranges to Staffli object 
## After filtering the dimensions of the experiment is: [27127 genes, 46278 spots]

Add a new column with the combined “protocol” and “ID”

SI$protocol_array <- paste0(SI$protocol, " : ", SI$arrayID)
ST.FeaturePlot(SI, features = "nFeature_RNA", ncol = 3, label.by = "protocol_array", show.sb = FALSE, pt.size = 1.5)

## Loading required namespace: viridis

Load H&E images

SI <- LoadImages(SI, time.resolve = FALSE, xdim = 1e3)

## Loading images for 14 samples: 
##   Reading ../data/spaceranger_output/smallintestine/V19T26-028_A1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 1 image from 1979x2000 pixels to 1000x1011 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V19T26-028_B1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 2 image from 1979x2000 pixels to 1000x1011 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V19T26-028_C1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 3 image from 1979x2000 pixels to 1000x1011 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V19T26-028_D1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 4 image from 1979x2000 pixels to 1000x1011 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-048_A1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 5 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-048_B1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 6 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-048_C1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 7 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-048_D1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 8 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-050_A1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 9 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-050_B1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 10 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-050_C1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 11 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10F24-050_D1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 12 image from 2000x1937 pixels to 1000x968 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V10S29-108_A1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 13 image from 1979x2000 pixels to 1000x1011 pixels 
##   Reading ../data/spaceranger_output/smallintestine/V11B18-363_D1/spatial/tissue_hires_image.png for sample 1 ... 
##   Scaling down sample 14 image from 1979x2000 pixels to 1000x1011 pixels

Apply rigid transformations to obtain a rough alignment of the tissue sections.

# Warp transform
SI <- WarpImages(SI, verbose = TRUE, transforms = list("3" = list(angle = -90), "4" = list(angle = -90), 
                                                         "5" = list(angle = 160), "6" = list(angle = 160),
                                                         "7" = list(angle = 20), "8" = list(angle = 20),
                                                         "9" = list(angle = 20), "10" = list(angle = 20),
                                                         "11" = list(angle = 20), "12" = list(angle = 20),
                                                         "13" = list(angle = -140), "14" = list(angle = 30)))

## Creating dummy masks ...Loading masked image for sample 3 ... 
## Warping pixel coordinates for 3 ... 
## Warping image for 3 ... 
## Warping image mask for 3 ... 
## Finished alignment for sample3 
## 
## Loading masked image for sample 4 ... 
## Warping pixel coordinates for 4 ... 
## Warping image for 4 ... 
## Warping image mask for 4 ... 
## Finished alignment for sample4 
## 
## Loading masked image for sample 5 ... 
## Warping pixel coordinates for 5 ... 
## Warping image for 5 ... 
## Warping image mask for 5 ... 
## Finished alignment for sample5 
## 
## Loading masked image for sample 6 ... 
## Warping pixel coordinates for 6 ... 
## Warping image for 6 ... 
## Warping image mask for 6 ... 
## Finished alignment for sample6 
## 
## Loading masked image for sample 7 ... 
## Warping pixel coordinates for 7 ... 
## Warping image for 7 ... 
## Warping image mask for 7 ... 
## Finished alignment for sample7 
## 
## Loading masked image for sample 8 ... 
## Warping pixel coordinates for 8 ... 
## Warping image for 8 ... 
## Warping image mask for 8 ... 
## Finished alignment for sample8 
## 
## Loading masked image for sample 9 ... 
## Warping pixel coordinates for 9 ... 
## Warping image for 9 ... 
## Warping image mask for 9 ... 
## Finished alignment for sample9 
## 
## Loading masked image for sample 10 ... 
## Warping pixel coordinates for 10 ... 
## Warping image for 10 ... 
## Warping image mask for 10 ... 
## Finished alignment for sample10 
## 
## Loading masked image for sample 11 ... 
## Warping pixel coordinates for 11 ... 
## Warping image for 11 ... 
## Warping image mask for 11 ... 
## Finished alignment for sample11 
## 
## Loading masked image for sample 12 ... 
## Warping pixel coordinates for 12 ... 
## Warping image for 12 ... 
## Warping image mask for 12 ... 
## Finished alignment for sample12 
## 
## Loading masked image for sample 13 ... 
## Warping pixel coordinates for 13 ... 
## Warping image for 13 ... 
## Warping image mask for 13 ... 
## Finished alignment for sample13 
## 
## Loading masked image for sample 14 ... 
## Warping pixel coordinates for 14 ... 
## Warping image for 14 ... 
## Warping image mask for 14 ... 
## Finished alignment for sample14

Add a new metadata column with a combined protocol and lung ID label

SI$paper_id <- "SI"

SI$protocol_id <- gsub(pattern = " ", replacement = "_", paste0(SI$protocol, "_ID: ", SI$paper_id))
ST.FeaturePlot(SI, features = "nFeature_RNA", ncol = 3, label.by = "protocol_id")

## Warning: Removed 174 rows containing missing values (geom_point).

Add manual annotations

SI <- ManualAnnotation(SI, type = "raw")

meta.data <- readRDS("../data/R_objects/SI_metadata_selections")
meta.data <- meta.data[colnames(SI), ]
SI@meta.data$labels <- meta.data$labels
SI@meta.data$labels[is.na(SI@meta.data$labels)] <- "background"

submucosa_13 <- rownames(subset(SI@meta.data, arrayID == "V10S29-108_A1" & labels == "background"))
SI@meta.data[submucosa_13, "labels"] <- "submucosa"

Filter out background

---

Keep only labeled spots.

SI <- SubsetSTData(SI, expression = labels %in% c("mucosa", "submucosa", "muscularis", "TLS", "serosa"))

Biotype content

---

Load hgenes.tsv containing biotype annotations and add the “mt_protein_coding” and “rb_protein_coding” categories.

ensids <- read.table("../data/genes/hgenes.tsv", header = T)
rownames(ensids) <- ensids$gene_name

ensids$gene_type[grep(pattern = "^MT-", x = ensids$gene_name)] <- "mt_protein_coding"
ensids$gene_type[grep(pattern = "^RPL|^RPS", x = ensids$gene_name)] <- "rb_protein_coding"

Figure 4 c - Biotype content as pie charts

---

First, we add a new column to our meta data with approximate storage times (time after sample collection). Then we calculate gene attributes: gene name, gene counts and biotype (gene_type). From this gene attribute data.frame we can then calculate percentages of each biotype for each storage time. Note that all biotypes that are not “IG_C_gene”, “IG_J_gene”, “IG_V_gene”, “TR_C_gene”, “TR_J_gene”, “TR_V_gene”, “lncRNA”, “protein_coding”, “rb_protein_coding” or “mt_protein_coding” are pur in a category called “other”.

# Add approximate storage time
SI$time <- SI[[]] %>%
  mutate(time = factor(case_when(seq_date == "200323" ~ "~ 1 month",
                          seq_date == "200923" ~ "~ 6 months",
                          seq_date %in% c("220404", "220506") ~ "~ 2 years"), 
                       levels = c("~ 1 month", "~ 6 months", "~ 2 years"))) %>% # ifelse(SI$seq_date %in% c("220404", "220506"), "220404-220506", SI$seq_date)
  pull(time)

# calculate gene attributes
gene_attr <- do.call(rbind, lapply(levels(SI$time), function(i) {
  x <- rowSums(GetAssayData(SI, slot = "counts", assay = "RNA")[, SI$time == i])
  data.frame(gene = names(x), count = x, gene_type = ensids[names(x), "gene_type"], time = i)
}))

# Calculate percentages for each storage time
gene_attr <- gene_attr %>% group_by(time, gene_type) %>% 
  summarize(N = sum(count)) %>% 
  group_by(time) %>% 
  mutate(p = N/sum(N)) %>%
  mutate(gene_type = factor(ifelse(gene_type %in% c("IG_C_gene", "IG_J_gene", "IG_V_gene",
                          "TR_C_gene", "TR_J_gene", "TR_V_gene",
                          "lncRNA", "protein_coding", "rb_protein_coding",
                          "mt_protein_coding"), gene_type, "other"), levels = 
                            c("protein_coding", "rb_protein_coding",
                          "mt_protein_coding", "lncRNA", "IG_C_gene", "IG_J_gene", "IG_V_gene",
                          "TR_C_gene", "TR_J_gene", "TR_V_gene", "other")))

## `summarise()` has grouped output by 'time'. You can override using the
## `.groups` argument.

gene_attr$time <- factor(gene_attr$time, levels = c("~ 1 month", "~ 6 months", "~ 2 years"))

# Create a color palette
cols <- c("#332288", "#6699CC", "#88CCEE", "#44AA99", "#117733", "#999933", "#DDCC77", "#661100", "#CC6677", "#882255", "#AA4499")

# Create pie charts
p <- ggplot(gene_attr, aes(x = "", y = p, fill = gene_type)) +
  geom_bar(width = 1, stat = "identity") +
  coord_polar("y", start = 0) +
  facet_wrap(~ time, ncol = 2) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.ticks = element_blank(), axis.text.x = element_blank(),
        strip.text = element_text(size = 14, colour = "black"), 
        legend.text = element_text(size = 14, colour = "black"), 
        legend.position = c(0.8, 0.3)) +
  labs(x = "", y = "", fill = "") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = cols)
p

pdf(file = "plots/biotype_content.pdf", width = 6, height = 6)
print(p)
dev.off()

We can plot the labels for the five regions as a spatial map.

ST.FeaturePlot(SI, features = "labels", ncol = 4, label.by = "time", show.sb = FALSE)

Supplementary figure - small intestine labels (five manually annotated regions)

---

# Create H&E image plots
plots_HE <- lapply(paste0(1:14), function(i) {
  im <- GetStaffli(SI)@rasterlists$processed[[i]]
  g <- grid::rasterGrob(im, width = unit(1, "npc"), height = unit(1, "npc"), interpolate = TRUE)
  ggplot() +
    annotation_custom(g, -Inf, Inf, -Inf, Inf) +
    theme_void() +
    theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0))
})

# Plot labels
labels_plot <- lapply(1:14, function(i) {
  p <- ST.FeaturePlot(SI, features = "labels", indices = i, show.sb = FALSE, pt.size = 0.7,
                 cols = c("mucosa" = "#AA4499", "submucosa" = "#DDCC77", "muscularis" = "#CC6677", "TLS" = "#77AADD", "serosa" = "#771155")) & 
    theme(plot.title = element_blank(), legend.position = "none", strip.text = element_blank())
  p <- ggrastr::rasterize(p, layers = "Point", dpi = 300)
})

# Combine H&E with label plots for each tissue section
plots <- lapply(1:14, function(i) {
  p <- plots_HE[[i]] | labels_plot[[i]]
  return(p)
})

# Patch together tissue sections
p1 <- cowplot::plot_grid(plotlist = plots[c(1:4, 13:14)], ncol = 1)

## Warning: Removed 2 rows containing missing values (geom_point).

p2 <- cowplot::plot_grid(plotlist = c(plots[5:8], list(NULL, NULL)), ncol = 1)

## Warning: Removed 77 rows containing missing values (geom_point).

## Warning: Removed 95 rows containing missing values (geom_point).

p3 <- cowplot::plot_grid(plotlist = c(plots[9:12], list(NULL, NULL)), ncol = 1)

# Create final patchwork
p <- p1 + p2 + p3 + patchwork::plot_layout(ncol = 3)
p

pdf(file = "../Suppl_figures/Suppl_Figure_smallintestine_labels/suppl_figure_HE_and_labels.pdf", width = 10, height = 10)
print(p)
dev.off()

Figure 4a - H&E image, top panel

---

Here we export the H&E image for the second tissue section. We also crop the image to remove the Visium fiducials.

par(mar = c(0,0,0,0))
plot(SI@tools$Staffli@rasterlists$raw$`2`[116:883, 95:860])

png(filename = "plots/HE.png", width = 860 - 95, height = 883 - 116)
par(mar = c(0,0,0,0))
plot(SI@tools$Staffli@rasterlists$raw$`2`[116:883, 95:860])
dev.off()

Figure 4a - labels as a spatial map (bottom panel)

---

And the manual annotations for the bottom half of Fig. 4a.

gg <- cbind(SI[[]], GetStaffli(SI)@meta.data)
gg <- subset(gg, sample == "2")
dims <- GetStaffli(SI)@dims[["2"]]

limits_list <- list("2" = c(x_start = 95*2, y_start = 116*2, x_end = 860*2, y_end = 883*2))

p <- ggplot(gg, aes(warped_x, dims$height - warped_y, color = labels)) +
  geom_point(size = 1.7) +
  theme_void() +
  scale_color_manual(values = c("mucosa" = "#AA4499", "submucosa" = "#DDCC77", "muscularis" = "#CC6677", "TLS" = "#77AADD", "serosa" = "#771155")) +
  scale_x_continuous(expand = c(0, 0), limits = c(limits_list[["2"]]["x_start"], 
                                limits_list[["2"]]["x_end"])) +
  scale_y_continuous(expand = c(0, 0), limits = c(dims$height - limits_list[["2"]]["y_end"], 
                                dims$height - limits_list[["2"]]["y_start"])) +
  theme(legend.position = "none")

p <- ggrastr::rasterize(p, layers = "Point", dpi = 300)
p

png(filename = "plots/annotations.png", width = (860 - 95)*2, height = (883 - 116)*2, res = 300)
print(p)
dev.off()

Figure 4b - distribution of unique genes across regions and time points as violin plots

---

p <- ggplot() + 
  geom_violin(data = SI[[]], 
              aes(arrayID, nFeature_RNA, fill = protocol), scale = "width") + 
  geom_hline(data = SI[[]] %>% group_by(protocol, time, labels) %>%
               summarize(Mean = round(mean(nFeature_RNA))), 
             aes(yintercept = Mean), linetype = "longdash") +
  geom_label(data = SI[[]] %>% group_by(protocol, time, labels) %>%
               summarize(Mean = round(mean(nFeature_RNA))), 
             aes("A", Mean, label = Mean), size = 5) +
  facet_grid(labels ~ time, scales = "free", space = "free") +
  scale_x_discrete(labels = c("", "Rep1", "Rep2", "Rep3", "Rep4", "Rep5", "Rep6", "Rep7", "Rep8")) +
  scale_y_log10() +
  scale_fill_brewer(palette = "Pastel2") +
  labs(x = "", y = "") +
  theme(panel.background = element_rect(fill = "white"), panel.grid = element_line(colour = "lightgray", size = 0.5), 
        strip.text = element_text(size = 14, colour = "black"), 
        axis.text = element_text(size = 14, colour = "black"), 
        legend.text = element_text(size = 14, colour = "black"),
        legend.title = element_text(size = 16, colour = "black"), 
        legend.position = "bottom")

## `summarise()` has grouped output by 'protocol', 'time'. You can override using
## the `.groups` argument.
## `summarise()` has grouped output by 'protocol', 'time'. You can override using
## the `.groups` argument.

p

pdf(file = "plots/unique_genes_overview.pdf", width = 14, height = 10)
print(p)
dev.off()

Normalize data

---

SISplit <- lapply(levels(SI$time), function(i) {
  cat(sprintf("Started processing data set %s", i), "\n")
  SIs <- SubsetSTData(SI, expression = nFeature_RNA > 100)
  SIs <- SubsetSTData(SIs, expression = time %in% i)
  SIs <- SIs %>% 
    NormalizeData()
  SIs <- SetIdent(SIs, value = "labels")
  return(SIs)
})

## Started processing data set ~ 1 month 
## Started processing data set ~ 6 months 
## Started processing data set ~ 2 years

DE test

---

Here we run the test for each storage time separately, while keeping only 1000 spots per group to speed up computations. We keep genes with a positive avg_log2FC value and with adjusted p-values lower than 0.01. Note that we test our manually annotated “mucosa” spots against the background spots.

de.markers.list <- lapply(SISplit, function(se) {
  de.markers <- FindMarkers(se, only.pos = TRUE, max.cells.per.ident = 1000, ident.1 = "mucosa") %>% 
    filter(p_val_adj < 0.01)
  de.markers$gene <- rownames(de.markers)
  de.markers$time <- unique(se$time)
  return(de.markers)
})

Now we can create the UpSet plot showing the numbers of detected DE genes for each storage time.

First, we combine all de markers from our de.markers.list list into one data.frame. Then we group the data.frame by gene and check whether the gene is present at one or more storage time points. Thereafter we calculate the number of genes detected in each category, where the categories are: “~ 1 month”, “~ 6 months”, “~ 2 years” or any combination of these time points. The expressionInput is a named integer vector where each element is named by the category described above.

de.markers <- do.call(rbind, de.markers.list)

expressionInput <- de.markers %>% 
  group_by(gene) %>%
  summarize(group = paste0(time, collapse = "&")) %>%
  group_by(group) %>%
  summarize(N = n()) %>%
  arrange(-N)
expressionInput <-  setNames(expressionInput$N, nm =  expressionInput$group)

p <- upset(fromExpression(expressionInput), order.by = "freq")
p

pdf(file = "plots/upset_DE.pdf", width = 6, height = 3)
print(p)
dev.off()

Figure 4e - spatial maps of selected enterocyte markers

---

The enterocyte markers were selected from the Human Gut Cell Atlas.

# Normalize the SI Seurat object
SI <- NormalizeData(SI)

# prepare data for plots
enterocyte.markers <- c("ANPEP", "RBP2", "DGAT1", "FABP2", "APOB")
gg <- cbind(SI[[]], GetStaffli(SI)@meta.data, FetchData(object = SI, vars = enterocyte.markers)) %>%
  tidyr::gather(variable, value, enterocyte.markers) %>%
  select(arrayID, time, warped_x, warped_y, sample, variable, value)

## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(enterocyte.markers)` instead of `enterocyte.markers` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

# Create spatial maps
plots <- lapply(enterocyte.markers, function(i) {
  subset(gg, variable == i) %>%
  subset(sample %in% c("2", "5", "13")) %>%
    ggplot(aes(warped_x, warped_y, color = value)) +
      geom_point(size = 0.3) +
      theme_void() +
      facet_grid(time ~ .) +
      scale_x_continuous(expand = c(0, 0)) +
      scale_y_reverse(expand = c(0, 0)) +
      scale_color_gradientn(colours = viridis::magma(n = 11, direction = -1)) +
      labs(title = i) +
      theme(legend.position = "top", 
            strip.text = element_blank()) -> p
  ggrastr::rasterize(p, layers = "Point", dpi = 300)
})

p <- cowplot::plot_grid(plotlist = plots, nrow = 1)
p

pdf(file = "plots/markers_spatial_enterocytes.pdf", width = 12, height = 8)
print(p)
dev.off()

## quartz_off_screen 
##                 2

Export H&E images

for (i in c("2", "5", "13")) {
  im <- SI@tools$Staffli@rasterlists$processed[[i]]
  png(filename = paste0("HE_", i, ".png"), width = ncol(im), height = nrow(im), res = 300)
  par(mar = c(0,0,0,0))
  plot(im)
  dev.off()
}

de.markers %>% mutate(Freq = as.numeric(table(gene)[gene])) %>%
  filter(Freq == 3) %>%
  group_by(gene) %>%
  summarize(Mean = mean(avg_log2FC)) %>%
  arrange(-Mean)

date

---

date()

## [1] "Tue Aug 30 11:17:10 2022"

Session

---

devtools::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.1.3 (2022-03-10)
##  os       macOS Big Sur/Monterey 10.16
##  system   x86_64, darwin13.4.0
##  ui       unknown
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Europe/Stockholm
##  date     2022-08-30
##  pandoc   2.19.2 @ /Users/ludviglarsson/miniconda3/envs/R4.1/bin/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package         * version  date (UTC) lib source
##  abind             1.4-5    2016-07-21 [1] CRAN (R 4.1.3)
##  backports         1.4.1    2021-12-13 [1] CRAN (R 4.1.1)
##  beeswarm          0.4.0    2021-06-01 [1] CRAN (R 4.1.0)
##  bit               4.0.4    2020-08-04 [1] CRAN (R 4.1.3)
##  bit64             4.0.5    2020-08-30 [1] CRAN (R 4.1.3)
##  bmp               0.3      2017-09-11 [1] CRAN (R 4.1.3)
##  broom             1.0.0    2022-07-01 [1] CRAN (R 4.1.3)
##  bslib             0.4.0    2022-07-16 [1] CRAN (R 4.1.3)
##  cachem            1.0.6    2021-08-19 [1] CRAN (R 4.1.1)
##  Cairo             1.6-0    2022-07-05 [1] CRAN (R 4.1.3)
##  callr             3.7.2    2022-08-22 [1] CRAN (R 4.1.3)
##  car               3.1-0    2022-06-15 [1] CRAN (R 4.1.3)
##  carData           3.0-5    2022-01-06 [1] CRAN (R 4.1.3)
##  cli               3.3.0    2022-04-25 [1] CRAN (R 4.1.3)
##  cluster           2.1.4    2022-08-22 [1] CRAN (R 4.1.3)
##  codetools         0.2-18   2020-11-04 [1] CRAN (R 4.1.3)
##  colorspace        2.0-3    2022-02-21 [1] CRAN (R 4.1.2)
##  cowplot           1.1.1    2020-12-30 [1] CRAN (R 4.1.3)
##  crayon            1.5.1    2022-03-26 [1] CRAN (R 4.1.2)
##  data.table        1.14.2   2021-09-27 [1] CRAN (R 4.1.3)
##  deldir            1.0-6    2021-10-23 [1] CRAN (R 4.1.3)
##  devtools          2.4.4    2022-07-20 [1] CRAN (R 4.1.3)
##  digest            0.6.29   2021-12-01 [1] CRAN (R 4.1.1)
##  dplyr           * 1.0.9    2022-04-28 [1] CRAN (R 4.1.3)
##  ellipsis          0.3.2    2021-04-29 [1] CRAN (R 4.1.0)
##  evaluate          0.16     2022-08-09 [1] CRAN (R 4.1.3)
##  fansi             1.0.3    2022-03-24 [1] CRAN (R 4.1.2)
##  farver            2.1.1    2022-07-06 [1] CRAN (R 4.1.3)
##  fastmap           1.1.0    2021-01-25 [1] CRAN (R 4.1.0)
##  fitdistrplus      1.1-8    2022-03-10 [1] CRAN (R 4.1.3)
##  fs                1.5.2    2021-12-08 [1] CRAN (R 4.1.3)
##  future            1.27.0   2022-07-22 [1] CRAN (R 4.1.3)
##  future.apply      1.9.0    2022-04-25 [1] CRAN (R 4.1.3)
##  generics          0.1.3    2022-07-05 [1] CRAN (R 4.1.3)
##  ggbeeswarm        0.6.0    2017-08-07 [1] CRAN (R 4.1.0)
##  ggplot2         * 3.3.6    2022-05-03 [1] CRAN (R 4.1.3)
##  ggpubr          * 0.4.0    2020-06-27 [1] CRAN (R 4.1.3)
##  ggrastr           1.0.1    2021-12-08 [1] CRAN (R 4.1.1)
##  ggrepel           0.9.1    2021-01-15 [1] CRAN (R 4.1.3)
##  ggridges          0.5.3    2021-01-08 [1] CRAN (R 4.1.3)
##  ggsignif          0.6.3    2021-09-09 [1] CRAN (R 4.1.3)
##  globals           0.16.0   2022-08-05 [1] CRAN (R 4.1.3)
##  glue              1.6.2    2022-02-24 [1] CRAN (R 4.1.2)
##  goftest           1.2-3    2021-10-07 [1] CRAN (R 4.1.3)
##  gridExtra         2.3      2017-09-09 [1] CRAN (R 4.1.0)
##  gtable            0.3.0    2019-03-25 [1] CRAN (R 4.1.0)
##  hdf5r             1.3.5    2021-11-15 [1] CRAN (R 4.1.3)
##  highr             0.9      2021-04-16 [1] CRAN (R 4.1.0)
##  htmltools         0.5.3    2022-07-18 [1] CRAN (R 4.1.3)
##  htmlwidgets       1.5.4    2021-09-08 [1] CRAN (R 4.1.1)
##  httpuv            1.6.5    2022-01-05 [1] CRAN (R 4.1.2)
##  httr              1.4.4    2022-08-17 [1] CRAN (R 4.1.3)
##  ica               1.0-3    2022-07-08 [1] CRAN (R 4.1.3)
##  igraph            1.3.4    2022-07-19 [1] CRAN (R 4.1.3)
##  imager            0.42.13  2022-03-07 [1] CRAN (R 4.1.3)
##  irlba             2.3.5    2021-12-06 [1] CRAN (R 4.1.3)
##  jpeg              0.1-9    2021-07-24 [1] CRAN (R 4.1.3)
##  jquerylib         0.1.4    2021-04-26 [1] CRAN (R 4.1.0)
##  jsonlite          1.8.0    2022-02-22 [1] CRAN (R 4.1.2)
##  KernSmooth        2.23-20  2021-05-03 [1] CRAN (R 4.1.3)
##  knitr             1.39     2022-04-26 [1] CRAN (R 4.1.3)
##  labeling          0.4.2    2020-10-20 [1] CRAN (R 4.1.2)
##  later             1.2.0    2021-04-23 [1] CRAN (R 4.1.0)
##  lattice           0.20-45  2021-09-22 [1] CRAN (R 4.1.1)
##  lazyeval          0.2.2    2019-03-15 [1] CRAN (R 4.1.3)
##  leiden            0.4.2    2022-05-09 [1] CRAN (R 4.1.3)
##  lifecycle         1.0.1    2021-09-24 [1] CRAN (R 4.1.1)
##  limma             3.50.3   2022-04-07 [1] Bioconductor
##  listenv           0.8.0    2019-12-05 [1] CRAN (R 4.1.3)
##  lmtest            0.9-40   2022-03-21 [1] CRAN (R 4.1.3)
##  magick            2.7.3    2021-08-18 [1] CRAN (R 4.1.3)
##  magrittr        * 2.0.3    2022-03-30 [1] CRAN (R 4.1.3)
##  MASS              7.3-58.1 2022-08-03 [1] CRAN (R 4.1.3)
##  Matrix            1.4-1    2022-03-23 [1] CRAN (R 4.1.2)
##  matrixStats       0.62.0   2022-04-19 [1] CRAN (R 4.1.3)
##  memoise           2.0.1    2021-11-26 [1] CRAN (R 4.1.1)
##  mgcv              1.8-40   2022-03-29 [1] CRAN (R 4.1.3)
##  mime              0.12     2021-09-28 [1] CRAN (R 4.1.1)
##  miniUI            0.1.1.1  2018-05-18 [1] CRAN (R 4.1.0)
##  munsell           0.5.0    2018-06-12 [1] CRAN (R 4.1.2)
##  nlme              3.1-159  2022-08-09 [1] CRAN (R 4.1.3)
##  openxlsx          4.2.5    2021-12-14 [1] CRAN (R 4.1.3)
##  parallelly        1.32.1   2022-07-21 [1] CRAN (R 4.1.3)
##  patchwork         1.1.2    2022-08-19 [1] CRAN (R 4.1.3)
##  pbapply           1.5-0    2021-09-16 [1] CRAN (R 4.1.3)
##  pillar            1.8.1    2022-08-19 [1] CRAN (R 4.1.3)
##  pkgbuild          1.3.1    2021-12-20 [1] CRAN (R 4.1.2)
##  pkgconfig         2.0.3    2019-09-22 [1] CRAN (R 4.1.0)
##  pkgload           1.3.0    2022-06-27 [1] CRAN (R 4.1.3)
##  plotly            4.10.0   2021-10-09 [1] CRAN (R 4.1.3)
##  plyr              1.8.7    2022-03-24 [1] CRAN (R 4.1.3)
##  png               0.1-7    2013-12-03 [1] CRAN (R 4.1.0)
##  polyclip          1.10-0   2019-03-14 [1] CRAN (R 4.1.3)
##  prettyunits       1.1.1    2020-01-24 [1] CRAN (R 4.1.0)
##  processx          3.7.0    2022-07-07 [1] CRAN (R 4.1.3)
##  profvis           0.3.7    2020-11-02 [1] CRAN (R 4.1.0)
##  progressr         0.10.1   2022-06-03 [1] CRAN (R 4.1.3)
##  promises          1.2.0.1  2021-02-11 [1] CRAN (R 4.1.0)
##  ps                1.7.1    2022-06-18 [1] CRAN (R 4.1.3)
##  purrr             0.3.4    2020-04-17 [1] CRAN (R 4.1.0)
##  R6                2.5.1    2021-08-19 [1] CRAN (R 4.1.1)
##  RANN              2.6.1    2019-01-08 [1] CRAN (R 4.1.3)
##  raster            3.5-21   2022-06-27 [1] CRAN (R 4.1.3)
##  RColorBrewer      1.1-3    2022-04-03 [1] CRAN (R 4.1.3)
##  Rcpp              1.0.9    2022-07-08 [1] CRAN (R 4.1.3)
##  RcppAnnoy         0.0.19   2021-07-30 [1] CRAN (R 4.1.3)
##  readbitmap        0.1.5    2018-06-27 [1] CRAN (R 4.1.3)
##  remotes           2.4.2    2021-11-30 [1] CRAN (R 4.1.1)
##  reshape2          1.4.4    2020-04-09 [1] CRAN (R 4.1.3)
##  reticulate        1.25     2022-05-11 [1] CRAN (R 4.1.3)
##  rgeos             0.5-9    2021-12-15 [1] CRAN (R 4.1.3)
##  rlang             1.0.4    2022-07-12 [1] CRAN (R 4.1.3)
##  rmarkdown         2.15     2022-08-16 [1] CRAN (R 4.1.3)
##  ROCR              1.0-11   2020-05-02 [1] CRAN (R 4.1.3)
##  rpart             4.1.16   2022-01-24 [1] CRAN (R 4.1.3)
##  rstatix           0.7.0    2021-02-13 [1] CRAN (R 4.1.3)
##  rstudioapi        0.14     2022-08-22 [1] CRAN (R 4.1.3)
##  Rtsne             0.16     2022-04-17 [1] CRAN (R 4.1.3)
##  sass              0.4.1    2022-03-23 [1] CRAN (R 4.1.2)
##  scales            1.2.1    2022-08-20 [1] CRAN (R 4.1.3)
##  scattermore       0.8      2022-02-14 [1] CRAN (R 4.1.3)
##  sctransform       0.3.4    2022-08-20 [1] CRAN (R 4.1.3)
##  sessioninfo       1.2.2    2021-12-06 [1] CRAN (R 4.1.1)
##  Seurat          * 4.1.1    2022-05-02 [1] CRAN (R 4.1.3)
##  SeuratObject    * 4.1.0    2022-05-01 [1] CRAN (R 4.1.3)
##  shiny             1.7.2    2022-07-19 [1] CRAN (R 4.1.3)
##  shinyjs           2.1.0    2021-12-23 [1] CRAN (R 4.1.3)
##  sp              * 1.5-0    2022-06-05 [1] CRAN (R 4.1.3)
##  spatstat.core     2.4-4    2022-05-18 [1] CRAN (R 4.1.3)
##  spatstat.data     2.2-0    2022-04-18 [1] CRAN (R 4.1.3)
##  spatstat.geom     2.4-0    2022-03-29 [1] CRAN (R 4.1.3)
##  spatstat.random   2.2-0    2022-03-30 [1] CRAN (R 4.1.3)
##  spatstat.sparse   2.1-1    2022-04-18 [1] CRAN (R 4.1.3)
##  spatstat.utils    2.3-1    2022-05-06 [1] CRAN (R 4.1.3)
##  stringi           1.7.8    2022-07-11 [1] CRAN (R 4.1.3)
##  stringr           1.4.1    2022-08-20 [1] CRAN (R 4.1.3)
##  STutility       * 0.1.0    2022-08-23 [1] local
##  survival          3.4-0    2022-08-09 [1] CRAN (R 4.1.3)
##  tensor            1.5      2012-05-05 [1] CRAN (R 4.1.3)
##  terra             1.5-21   2022-02-17 [1] CRAN (R 4.1.3)
##  tibble            3.1.8    2022-07-22 [1] CRAN (R 4.1.3)
##  tidyr             1.2.0    2022-02-01 [1] CRAN (R 4.1.3)
##  tidyselect        1.1.2    2022-02-21 [1] CRAN (R 4.1.3)
##  tiff              0.1-11   2022-01-31 [1] CRAN (R 4.1.3)
##  UpSetR          * 1.4.0    2019-05-22 [1] CRAN (R 4.1.3)
##  urlchecker        1.0.1    2021-11-30 [1] CRAN (R 4.1.3)
##  usethis           2.1.6    2022-05-25 [1] CRAN (R 4.1.3)
##  utf8              1.2.2    2021-07-24 [1] CRAN (R 4.1.0)
##  uwot              0.1.14   2022-08-22 [1] CRAN (R 4.1.3)
##  vctrs             0.4.1    2022-04-13 [1] CRAN (R 4.1.3)
##  vipor             0.4.5    2017-03-22 [1] CRAN (R 4.1.0)
##  viridis           0.6.2    2021-10-13 [1] CRAN (R 4.1.1)
##  viridisLite       0.4.1    2022-08-22 [1] CRAN (R 4.1.3)
##  withr             2.5.0    2022-03-03 [1] CRAN (R 4.1.2)
##  xfun              0.32     2022-08-10 [1] CRAN (R 4.1.3)
##  xtable            1.8-4    2019-04-21 [1] CRAN (R 4.1.0)
##  yaml              2.3.5    2022-02-21 [1] CRAN (R 4.1.2)
##  zeallot           0.1.0    2018-01-28 [1] CRAN (R 4.1.3)
##  zip               2.2.0    2021-05-31 [1] CRAN (R 4.1.0)
##  zoo               1.8-10   2022-04-15 [1] CRAN (R 4.1.3)
## 
##  [1] /Users/ludviglarsson/miniconda3/envs/R4.1/lib/R/library
## 
## ──────────────────────────────────────────────────────────────────────────────